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Chapter 1 


INTRODUCTION 

In the operation of a liquid-propellant system the injected 
propellants are converted by many physical and chonical processes into 
hot burned gases vMch are accelerated to si^rscmc velocity by passing 
through a converging- diverging nozzle. Since the operation of such a 
systan, however, is seldom perfectly smooth, the oscillations can be 
of either a destructive or nondestructive nature. Iftidestructive 
isisteadiness is characterized by randcmi fluctuations in tl« flow 
properties and includes the phenomena of turbulence and coni>iistion 
noise. Unsteady operation of a destructive nature, on the other hand, 
is characterized by organized oscillations in \diich there is a definite 
correlation between the fluctuations at two different locations in the 
conbustor. Such oscillations have a definite frequency and result in 
additional thermal and mechanical loads that the syst«n must withstand. 

In liquid-propellant systems experiencing unstable canbustion, 
heat transfer rates to the walls considerably exceed the corresp<mding 
steady state heat transfer rates, resulting in bum-out of the walls. 

If the system can survive these effects, mechanical vibrations in the 
system can cause mechanical failure or destroy the effectiveness of 
the delicate control and guidance systems. 

The phenanenon of cotrtjustion instability depends strongly iqjon 
the unsteady behavior of the confcustion process. The organized 
oscillaticHis of the gas within the system must be coupled with the 
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combustion process in such a my as to font a fee^ack loop. An 
understanding of this ccaipling between the coRd>ustion process and the 
wave motion is necessary in order to predict the stability characteristics 
of the systan. 

There are three cases of oscillations in the combustion 
instability problems. The most in^rtant form of combustion instability 
is known as high frequency instability. As the name suggests, this 
type of instability represents the case of forced oscillations of the 
conbustion chamber gases which are driven by t}% unsteady combustion 
and intract with the rescmance properties of the conibustion geanetry. 

The observed frequencies, vhich are as high as 10,000 cycles per second, 
are very close to those of the natural acoustic modes of a closed-ended 
chamber of the same geometry as the one experiencing unstable cond)ustion. 

High frequency canbustion instability can resemble any of the 
following acoustic modes: (1) longitudinal, (2) transverse, and 

(3) combined logitudinal- transverse modes. Longitudinal oscillations 
are usually observed in chambers whose length to diameter ratio is 
much greater than one; in this case the velocity fluctuaticms are 
parallel to the axis of the chamber and the disturbances depend only on 
one space dimension. For much shorter chanhers the transverse mode of 
instability is most frequently observed. Transverse oscillations in 
rocket motors are characterized by a ccmiponent of the velocity' 
perturbation which is perpendicular to the axis of the chamber and a 
three-dimensional disturbance field. Such oscillations can take either 
of two forms: (1) the standing form in which the nodal surfaces are 
stationary, and (2) the spining form in which the nodal surfaces rotate 
in either the clockwise or ccwnterclockwise direction. 
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Hm theoretical study of iKsilinear vibration ai»i conbustion 
instability appears to have beoi initiated by a British group [1] in 
1940. At that tiae the group testing a small solid-pio^iellant rocket 
motor (^served sudden increases of pressure to twice the e^qpected 
level, enough to destroy a motor of flight weight. Since the early 
1950 's nuch experimental and analytical research has been developed to 
better understand the pheiuneiKm of high frequency conbustion 
instability. Most of the theories presented prior to 1966 were 
restricted to circimistances in \diich the amplitudes of the pressure 
oscillations were infinitesimally small, that is, in tl^ linear regime. 
The case of longitudinal instability was studied by Crocco [2] as well 
as the studies of transverse instability by Scala [3] , Reardon [4] , 
Culick [5], and Zinn [6]. 

In the field of finite an^litude (ntxilinear) ccxnbustion 
instability, mathematical difficulties have precluded any exact 
solutions, and approximate methods and numerical analyses have been used 
almost exclusively. For this reason publications in this field are 
relatively scarce. Notable among these is the work of Maslen and Moore 
[7] vdio studied the behavior of finite anplitude transverse waves in a 
circular cylinder. Their major conclusion was that, unlike logitudinal 
oscillations, transverse waves do not steepen to form shock waves. 

Maslen and Moore, how-ewr, considered only fluid medvinical effects; 
they did not consider the influences of the combustion process. 

One of the first nonlinear analyses to include the effects of 
the combustiCTi process and resulting steady state flow was performed bv 
Priem and Guentert [8]. In this investigatiai the problem was made 
one- dimens i<mal by considering the behavior of tangential waves 
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traveling in an annular section of the conbustor of a liquid prc^llant 
rocket motor. In more recent years other investigators such as 
Burstein [9] have attespted to solve manerically the equations 
describing instabilities that depend two space dimeisions. Although 
the resulting solutions resonble experimentally (^served combustion 
instability, this method requires excessive conputer time, aiKl studies 
of this type for three-dimensioial oscillations will have to await the 
develqnnent of a much faster breed of oxqmters. Ptxiiell [10] studied 
nonlinear coirbustion instability in liquid propellant rocket engines 
using the Galerkin method. He concluded that the following ncmlinear 
mechanisms were inportant in determining the nonlinear stability 
characteristics of the syst«n; (1) the transfer of energy between modes, 
(2) the self-ccHipling of a mode with itself, and (3) a nonlinear 
COTibustion mass source. Also Powell's work indicates that the Galerkin 
method can be successfully applied in the soluti<m of mnlinear 
cranbustion instability problans. The numerical solution of the 
resulting ordinary differential equa*^ions will, in most cases, require 
much less computer time than a finite difference soluticm of the 
original partial differential equations. 

The present work is motivated by a desire to determine whether 
or not srane of the features of Powell's numerical solutions can be 
reproduced analytically by means of perturbation methods. For 
simplicity a thin rectangular chamber and a thin annular chanber are 
ccMisidered rather than the full cylinder investigated by Powell [10]. 
This eliminates wie space coordinate but still allows sufficient 
generality to make possible the investigation of both mode-coupled 
annular and self-coupled rectangle instabilities. Also for simplicity 


a siinplified equation of state for the gas and a sinplified description 
of the interphase mass transfer are enployed. Tlwse simplifications 
make it unnecessary to explicitly consider t)^ gas-phase energy equation 
and the fuel-drop phase balance laws. The unsteady condxistion response 
is characterized by a simple time-delay mechanisa. 

The Galexidn method is used to solve tlm partial differential 
equation governing the velocity potential. This produces a set of 
nonlinear ordinary differential equations which are solved by singular 
perturbation metlnds. This alloms limit cycles to be determined 
analytically and the approad^s to the limit cycle to be found from 
sinple numerical integrations. 


CSiapter 2 


GOMBUSTIGN EQUATIONS 

It is the objective of tl% discussions presented in this 
chapter to analyze the ccndsustor governing equati<»is and redtice than to 
an analytical form that is amenable to mathematical solution. The 
sinplification will be done in a manner that will allow thr resultant 
equations to preserve both the mathematical and physical essoice of the 
original problon. 

The ccmbustion chamber medim will be considered as a two-phase 
mixture of liquid propellant droplets and gaseous condjusticn products. 

The primary reaction between the fuel and oxidizer is assumed to occur 
inmediately iq>on vj^orization of the propellant drqjlets; that is, the 
vaporizaticm process is rate craitrolling. Therefore, the gas phase can 
be treated as a single constituent. For sijiplicity the equation of 
state is assuned to have the form of a linear relation between pressure 
and density. The presence of burning liquid droplets in the chamber will 
be represented as a continuous distribution of mass sources for the 
gas. The burning process will be treated as history dependent, this 
dependence being represented by a suiple time-delay mechanism. 

The following equaticMis are derived by applying the laws of 
conservation of mass, and momentum to an aibitrary stationary control 
volune. (For similar analyses see Readar [4], and Culick [5].) Because 
of the use of the sinplified equation of state, explicit consideration 
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of tlw energy ec)uatian is not required. ‘HMy are presmted first in 
dinmsional fbm aiKl %dll be non-dinensionalized later. 

Balance of mass: 

^ (pV) - B* (2.1) 

Balance of momoitun: 

p*(0t^ ♦ 0* • ^*0*) • -^V (2.2) 

where in equatxm (2.1), and (2.2) p , P , B , and 0* are tl« density, 
pressure, buming rate, and velocity in the chamber. Also v is the 
dunensioial gradient operator, and a coona followed by a s'jhscripted 
independent variable denotes the partial derivative with respect to that 
variable (in this case time) . 

The idealized equation of state takes on the mathematical form 

a A 

P -a^p (2.3) 

vdiere a* is the speed of sound at the reference state. Substituting 
(2.3) into (2.2), > elds 



(2.4) 


The governing equations will now be nondimens ionali zed with 
respect to a steady reference state, which will be denoted by the 
subscript "r." Usually this reference state is taken to he the 
stagnation condition at the injector face, but if the variation of 
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steady state pit^rties is iMgli{ible the referatce condition can be 

the corresponding constant steady state conditiwi. All lengths will be 

* 

referred to sane characteristic length L . (Sudi as the chanber radius 
of a cylindrical chasi)er.) The characteristic velocity is t)w speed of 

smnd at tlM reference state, and the cham:teristic time is the vnve 

* * 

travel time a^ /L . The dimensionless quantities may then be written 
as: 


-k* 

0 - 


a_ 


(2.5) 


( 2 . 6 ) 




t - 4 t 


(2.7) 

( 2 . 8 ) 


P “ -tT 




(2.9) 


8-6 


P a 

r 


( 2 . 10 ) 


Combining equaticwis (2.1), (2.3), and (2.4) with equations (2.S-2.10), 
the dimensionless governing equatiwis become 


p,^ ♦ ^ • (pi5) ■ 8 
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( 2 . 11 ) 

0,^ ♦ 3 • ^ ■ 4(log p) (2.12) 

p ■ p (2.13) 

>^re "log" neans natural logarithm. 

It will be considerable iMlp in the solution of the system of 
governing differential equations (2.11) and (2.12) if a first integral of 
t}M moraentim equation can be fovnd. This can be accomplished if the flow 
is irrotational, because the velocity vector can then be expressed as the 
gradient of a scalar potential. Uhder these conditions it is possible 
to combine the equations to obtain a single partial differential 
equation governing the behavior of tlw velocity potential. The pressure 
is then related to the velocity potential by the integrated raomenttra 
equation. To prove that flow is irrotational, it can be done by taking 
"curl" from equation (2.12) and also lising the vector identities 

^ • (^xF) • 0 

?x(^A) - 0 

F") - ? • ^ 

to find 


DW 

Dt 


- ^ - (^ • 3) i5 


(2.14) 


where W is vorticity vector, and defined by 
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^ M (2.15) 

Equation (2.14) can be satisfied by 

if - 0 (2.16) 

vdiich vrill be true provided 

W(r,0) - 0 (2.17) 

therefore equation (2.16) shows that flow is irrotational. 

Introducing the velocity potential defined by 

0 = ^ij) (2.18) 

substituting (2.18) into (2.12) and (2.11) 

-('t'-t * * 

p = e'- 

p,^ + + V4» • 

Due to their nonlinear natures the system of wave equations 

(2.19) and (2.20) cannot be solved exactly. In order to obtain simpler, 

but approximate, form, an order of magnitude analysis will be performed to 

detemine which terms can be neglected. In this approximation schane 

each perturbation quantity will be assigned the order of magnitude 0(e) 

where e is a measure of the magnitude of the initial disturbance. It 

then follows that products of perturbation quantities are of second 
2 

order or 0(e ). To obtain an equation that is correct to a given order 


yields 
p = B 


(2.19) 

( 2 . 20 ) 


u 


all terms of higher order are iwglected. By this definition it is se«i 
that first order equations are necessarily linear. 

The equations derived from this perturbation analysis are 
expected to be valid as long as the anq[)litude of the perturbation 
quantities is finite l»it smaller than unity. If in addition it is 
assumed that the Mach nuaber of the steady state flow is small, 
additicmal tenns can be neglected. The behavior of the steady’ state 
quantities is governed by the time independent versions of equatirais 
(2.19) and (2.20). Setting the time derivatives equal to zero, and 
retaining only terms of 0 (1) and 0(e) yields the steady-state equations. 
Sidjstituting into these 


'j' ■ el^(2) 
6 - F(2) 

f - 0(c) 


and letting 


3 ■ w ■ CO 


yield 


P ° I - 



( 2 . 21 ) 


( 2 . 22 ) 

(2.23) 


(2.24) 


(2.25) 



(2.26) 


For c<<l (2.26) reduce to 


A 
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d% - 

5? ' C2.27) 

It will be assuned that the cQmbusti(»i is udfomly distributed, w is 
ccwistant, therefore (2.27) can be integrated to find steady state 
velocity distribution 

(J ■ oz + c (2.28) 

U(0) - 0 (2.29) 

thus 

U ■ az (2.30) 

The ccnplete velocity potential is now defined by sum of steady 
state and unsteady state contribution, and a second order analysis will 
be performed on the simplified wave equation [i.e. eqs. (2.19-2.20)] to 
obtain a set of second order equation involving the velocity potential. 
After a first integral of the momentum equation is obtained these 
equations will be combined to yield a single nonlinear differential 
equation governing the behavior of the velocity potential. To this end 
let 


e(? + 4>(x,y,z,t)) 

(2.31) 

w + ew ■ e(o ♦ eo) 

(2.32) 


substituting (2.31) into (2.19) yields 
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p - 1 - - (o^ ♦ • ^0) - 2?0,J ♦ -— (2.33) 

substituting (2.33) and (2.31) into (2.20) yields 

♦ .tt ■ ♦ e[ 0 (|>,^ + 2oi|),^^ + 2 % • ^(♦f^) + o] " 0 (2.34) 


note that 

♦ “ v2«|> + 0(e) (2.35) 

It is to be noted that equation (2.34) has the form of an 
inhomogeneous wave equatioi. 

Substituting (2.35) into (2.34) yields 

4),^^ - * e(o<l>,t + 2 % • + a) « 0 (2.36) 

To roughly account for changes in the frequency spectrum due to 
baffles, cavities, etc., the governing equatioi will be modified as 
follows 


* Zh • + a) = 0 (2.37) 

where the term multiplied by k is a correction factor. This can also be 
written as 


4>,tt ■ ^*’4' + 20(^4^ • ^(4>»^) - *c^‘’(4>,^,^)) = -w 


(2.38) 


Chapter 3 


CTMBUSTION IN A THIN RECTANGULAR 
CHAMBER 

In this chapter the modified form of the governing equation 
(2.36) will be used to study the effect of nonlinearities on the 
stability of a thin rectangular chanber in which the distributed 

I 

combustion process is that described by Croco [2] . For simplicity one 
dimensional instabilities will be considered for a ccmbustor geometry 
with uniformly distributed injection (see Fig. 1). In this case all 
lengths are referred to characteristic length L*, the chamber length. 

The dimaisionless coordinates may then be written as: 

* 

X = ^ (3.1) 

L 

* 

y = ^ (3.2) 

L 

* 

(3.3) 

L 

Therefore, the nondimesional combustor geometry is given by 

Fig. 2. 

To formulate the combustion function it is first necessary to 
recall that the dimensionless pressure p>erturbation can be expressed as: 
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P - Pq 
Po® 
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♦ 0(e) 


(3.4) 


The oxnbustion process is assuaed to be described by t)w 
pressure- sensitive tine-delay model devel(^>ed by Croccc [2] and employed 
by Powell [10]. 

In the presoit notation this can be written as: 


o - an(q - qj.),q^ ■ q(t - T) (3.5) 

vhere n is called the interaction index (a measure of the stroigth of 
the con^ustion process) and T is called t}» (dimensionless) time-delay 
parameter (a measure of the ii^ortance of history effects) . 

The governing equation (2.36) now becomes: 


^»tt ■ * e(o((l - n)^,^ + n(<|>^),^) + 2^0 • ^(0,^.) 

-<7^(0, ^^)) + O(e^) - 0 (3.6) 

It is now necessary to select a potential function to predict 
the behavior of the amplitutes of the various acoustic modes of the 
chamber. As a first approximation attention will be restricted to the 
first term of a Fourier series, satisfying the boundar\' condition at 
y = ±7t/2, that is 

0(y»t) - F(t)siny (3.7) 

Substituting (3.7) into (3.6) (and letting k = 0 for sinplicity) 


yields 


16 

+ F + e(o(l - n)F + + 8/3irF?) + O(e^) - 0 (3.8) 


i^re 


(*) - d( )/dt 

Equation (3.8) has the same mathematical fbim as the equation 
obtained by Powell [10] in his study of purely radial instability in a 
cylindrical comlHistion chamber. His numerical solutions showed that 
the stability of the solution depended on the initial anplitude of the 
disturbance. Powell [10] also found some cases which were 
classified as stable by the linear theory were predicted to be unstable 
by the nonlinear theory and vice-versa. In the present work a 
combination of the two-variable and strained-parameters perturbation 
methods (see, for instance Nayfeh [11]) will be used to obtain 
approximate closed- form solution to (3.8) vrfiich will be used to explain 
the behavior discovered by Powell [10]. 

Let, 

■ 8e/3iT (3.9) 

To begin the two- variable perturbation process, let 

^ - t (3.10) 

n » et (3.11) 

and subtitute (3.9), (3.10), and (3.11) into (3.8) to get: 

♦ F ♦ e.FF,. + e(2F,. o((l - n)F,. ♦ n(FT>), )) + O(e^) » 0 
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Assuming 

F - Fq ♦ eFj + O(e^) (3.13) 

the governing equation for first and second order ai^roximations can be 
readily vrritten as: 




+ F„ 




(3.14) 


• + o((l - n)pQ,^ ♦ nCpQ^),^)) 

♦ O(e^) - 0 (3.15) 

An approximate solution to (3.14) will now be found by using the 
method of strained-parameters with e as the perturbation parameter. 
Toward this end let 


C - (1 + ♦ — -)2 

substituting (3.16) into (3.14), and assuming 

*^0 " ^00 ^ h^’oi ^ h\z * 

yields 

%’zz ^00 ” ® 
^Ol’zz * ^01 * '%%'z 


(3.16) 


(3.17) 


(3.18) 


(3.19) 
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^02'zz * ^02 “ ■^00*’01»z ■ *’oAo»z * ^l%»zz 


(3.20) 


It is assi;roed that tlw solution equation (3.18) can be written 


Fqq - C(n)cos(z - a(n)) 


(3.21) 


wheie C(n) , and ct(n) are the anplitude arei phase angle respectively. 

Substituting (3.21) into (3.19), and finding a particular 
solution of the resulting equation y-elds 


Fqi ■ -1/6 C sin2(z - a) 


(3.22) 


and substituting (3.21) and (3.22) into (3.20) yields 

Fq 2 ,^^ + Fq 2 ” (1/12 - 2WpCcos(z - a) + 1/4 Aos3(z - a) (3.23) 

The particular solution corresponding to the first term on the 
right hand side will be proportional to zsin(z - a). Thus limit 


as z goes to «> (violating the basic assunption of the perturbation 
process, namely that g^Fq^ and gJFq^ represent small correction to Fqq) 
unless this term vanishes. Thus 


Wj - 1/24 C" 


(3.24) 
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Substituting (3.24) tmd (3.16) into (3.21) yields 


2 ^ 2 . 


Fqq - C cos((l - 1/24 eJC‘)5 - a) + O(e^) 


(3.25) 


Note, 


pQ " ^'OO * °^®1^ 


(3.26) 


Novr by substituting (3.26) into (3.15) and using the fact that: 


sin(A - B - C) - sin(A - B)cosC - cos (A - B)sinC 


one finds that 


^1*5C * 


(2^ ♦ o((l - n)C + nOcosCl - ^V)T))- 


sinCd - ^V)C - a) - C(2^ + ohsin(l - j^eV)!)- 


cos((l - - Qt) 


(3.27) 


Eliminating singularity and using (3.9) and (3.11) yields 


C ♦ |€aC(l - n(l - cos(l 


8 2 

® “CnT)) =« 0 


27tt" 


(3.28) 


a + -^on sin(l - 2 €‘'C“)T = 0 (3.29) 

27tt 

The linear stability analysis can be recovered from (3.28) and 
(3.29) by neglecting the nonlinear terras in the arguments of the 
trigonometric functirais. This yields 
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C + JeaC(l - n(l - cosT)) - 0 


d + ^ohsinT ■ 0 


(3.30) 

(3.31) 


with initial conditions 


C(0) - 1 


o(0) » 0 


(3.32) 


(3.33) 


The appropriatf! closed- fora soluti<ms can be written as: 


C « ■ "Cl ■ cosT))t/2 


a • - 2 (ecmsinT)t 


It can be seen that perturbation in C will decay for 


(3.34) 

(3.35) 




(3.36) 


thus 


n » 


1 

1 - c»sT 


(3.37) 


is the neutral stability curve according to the linear theory (Fig. 3j. 

Equations (3.28) and (3.29) must be solved numerically. To get 
an idea of the behavior that will be predicted by numerical analysis, 
it is instructive to obtain approximate closed- form solutions to (3.28) 
and (3.29). This can be done by e>q)anding the trigonometrical functions 
appearing therein and retaining only the first teras to get 
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1 — 

C ♦ ^oC(l - n(l 


(cost ♦ — V^siiO'C^))) 
ii-r 


(3.38) 


o ♦ ^on(sinT - — ^jc^cosTC^) ■ 0 (3.39) 

ZTir* 

The initial conditions remain 

C(0) - 1 (3.40) 

a(0) - 0 (3.41) 

Eauations (3.38) aivt (3.39) may be written as the form 


C - aC ♦ bC^ • 0 (3.42) 

2 

a + eC^ + i - 0 (3.43) 

where a, b, e, and i are 

a - -^0(1 - n(l - COST)) (3.44) 

b » — l^^TsinT (3.45) 

llv 

e - -bcotT (3.46) 

i - ^ohsinT (3.47) 


The solutions of Equations (3.42) and (3.43) can be written 


as: 
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C(t) 


^at 

■ e 




b(l 




at^ 


(3.48) 


Inti - |(1 - - it (3.49) 

Information about the behavior of the solutions obtained in this 
duster is conveniently illustrated by considering stability bcaaidaries 
plotted in the n - T plane. By setting 


C - 0 


(3.50) 


in (3.28), the equaticm 


"e 


1_ 

1 - cos(l - 


8e^ 


27tt 



(3.51) 


is obtained. This is the ncHilinear neutral -stability curve. For pairs 
of values of n and T satisfying this equation the amplitude of the 
solution will remain fixed for all values of time. If e is equated to 
zero in (3.51) the linear neutrai-stability equatiwi 


1 

“ 1 - cost (3.52) 


gr 

is obtained. If (3 51) is expanded for «1 one obtains the 

27tt" 

approximate neutral-stability curve 


1 - cosT 


8c‘ 


27tt' 


TsinT 


(3.53) 
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This corresponds to the ai^roxijnate solution (3.48). Figures 4 
and S s\km plots of (3. SI) and (3. S3) respectively for the same values 
of e. It can be se«i that for c^l the two sets of curves are in close 
agreemmit. In fact, mraerical coqHitations have s)K>im that for 
equations (3. SI) and (3. S3) produce results that are practically 
indistinguishable in the region 0<T<2II. This suggests that the 
ai^roxinate solution (3.48) is likely to be a good approximation to the 
exact numerical results for moderate values of e. In \diat follows the 
validity of this hypothesis will be demonstrated. 

To an exaggerated scale typical curves plotted in Figure 4 looks 
as follows. 

n 3 


V 


1 


e 


The stability diagram can be conveniently divided into six 



regions as shown above. There are 
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I: 

T<it 

, n<nL 

II: 

T<ir 

, iY<n<n£ 

III: 

T<n 

» n>ng 

IV: 

T>tt 

, n<ng 

V: 

T>it 

» n£<n<nL 

VI: 

T>ir , 

, n>n^ 


The approximate solution (3.48) will now be used to infer the 
behavior of the solutirai in each region. Then some exact numerical 
results will be examined. The approximate equations are repeated below: 

at 

c = e (3.54) 

/l - |(1 - 

a = ^a(n(l - cosT) - 1) = - 1) (3.55) 


4g 

b = -^S_^TsinT 
27tt^ 


(3.56) 


The stability bomdary corresponds to a = b. 

In region I a<0, b>0. Substituting a = -a^ into (3.54) yields 


C 


■a^^t 


/ T 

1 * I (1 - e 1 ) 
^1 


(3.57) 


It can be seen that 
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lim C ■ 0 

(3.58) 

t -»■ oo 

Th\is the aaplitude decays to zero in region I. In region II 
a>0, b>0 and a<b. In region III a>0, b>0 and a>b. In either case 
(3.54) iraplles 


Ita c - /l" 

t -► 00 


(3.59) 


It can be seen that in both region II and region III a 
limit -cycle is sqjproached. For regicm II the anqplitude of the 
limit-cycle is less than that of the initial amplitude while for region 
III the limit-cycle anplitude is greater than the initial aiUJlitiKie. 

For a * b (3.54) reduces to 


C - 1 (3.60) 

which is the solution associated with the nonlinear stability curve. 

In region IV a<0, b<0 and a<b. In region V a<0, b<0 and a>b. 
Substituting a = and b * -b^^ into (3.54) me obtains 


C » 


e*V 



-2a,'t 
:p(l - e M 
‘^1 


(3.61) 


where for region IV a^^b^ and for region V a^<bj^. For the former region 
it is possible to take the limit as t goes to infinity to get 
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lin» C - 0 

(3.62) 

t ■+ 00 

In region IV, therefore, the amplitude decays to zero. For the 
latter region, <ai the other hand, the corresponding limit can not be 
taken because C becones infinite in a finite time. Denoting this time 
by t«. it can be seen fran (3.61) that 


b, -2a,t 

1 - ^(1 - e 1 ” ) - 0 or 




If ■ bj (3.63) reduces to 


(3.63) 


C - 1 


(3.64) 


v^rifiich corresponds to neutral stability. 

In region VI a>0 and b<0. Sidistituting b = -bj into (3.54) 

yields 


C 






3it 
(e - 1) 


( 3 . 65 ) 


In this case C bcconK's infinite when t = t^, where 


bi 2at,„ 

1 - — -(e - 1) =• 0 or 

vl 




1 

3a 


ln(l + 



(3.66) 


27 


Hear linear stability curve |a|«l. Expanding (3.63) and 
(3.66) in this region one finds that 





(3.67) 




(3.68) 


Thus, 


(t ) < (t ) . 

eo VI “V 


The results obtained in regions V and VI violate the ccmdition 
of validity of (3.54), namely that the magnitude of C be moderate. Thus 
the infoimation obtained abcHit the behavior of the solution in this 
region can be used only in a qualitative sense. Nevertheless, as vdll 
be seen, (3.61) and (3.65) correctly predict the main feature of exact 
solution. 

To sunmarize, the following behavior has been found based on the 
approximate solution (3.53). Below curve .4CE (regions I and IV) all 
disturbances decay to zero. Above cui-\'e ACF. the behavior of the 
solution depends on the value of T. For T<tt all distmionces evolve 
into limit cycles. In region II the value of the final ain^iUtude is 
less than that of the iiiitial amplitude, while in region III the 
situation is reversed. For T>it all disturbances increase without bound. 
The rate of increase is slower for region V than for i-egion VI. 
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Some typical nuaerical results obtained by solving equation 
(3.28) by a fourth*order Runge-kutta netl^ are shown in Figures 6-9. 
Figure 6 shows three sets of data c^tained along the line n - 1 in the 
vicinity of T ■ ^ (the curves are labeled in degrees rather than radians 
to better bring out the small differences in the values of T that were 
used). Each curve is represaitative of a different parametric region 
which included in the labeling of each curve. It can be seen that 
these results are in qualitative agreement with the predictions of the 
^proximate solution. Other results, not shown, reveal that the 
anplitude of the limit-cycle in region III increases rapidly as T 
increases. Such large limit-cycles are probably inconsistent with the 
order of magnitude assumption inherent in the perturbation method and, 
therefore, will not be discussed. 

Figure 7 shows three sets of data obtained along the line T ■ y 
in the vicinity of n = 1. The results for regions I and III are again 
in caiplete qualitative agreement with the jpproximate solution. It can 
be seen that the longest decay time is exhibited by data correspcHiding 
to a point on the borcteriine between regions I and II. The corresponding 
result based on the approximate analysis is found by taking the limit 
of (3.54) as a 0 to get 


C = 



(3.69) 


Since the decay in this case is algebraic rather than exponential (as it 
is for a 0) this case will exhibit the slowest decay in agreement with 
the exact solution. 
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Figure 8 shows three sets of data obtained along the line 
n ■ 1 in the vicinity of T • ^ . The case corresponding to region IV 
decays to zero in agreonent with the aj^roximate solution. In regions 
V and VI, however, the disturbances evolve into limit cycles rather 
than growing without bound as predicted by the approximate analysis. In 
contrast to the situation for T<tt, it does not appear possible to find 
any values of T for v^ch the limit-cycle anplitude will have a moderate 
value. Since the accuracy of solutions exhibiting large limit -cycle 
anplitude is doubtful, the difference between the unbcxmded growth 
preuicted by the approximate solution and the growth to a very large 
limit cycle predicted by the exact analysis is probably uninportant. It 
can be seen that the rate of growth of the solution increases as one 
proceeds from region IV to region VI as predicted by the approximate 
solution. 

Figure 9 shows two sets of data obtained along the line T » 
in the vicinity of n = 1. Again the differences in behavior between 
regions IV and VI are clearly illustrated. 

Powell [10] observed similar behavior to that discussed above 
in his numerical solutions to an equation similar to (3.28). The 
perturbation analysis given here makes it clear that the main effect 
causing the observed behavior is the frequenc- shift due to nonlinear 
amplitude -frequency dependence and the accompanying distortion of the 
neutral stability curves in the n - T plane. 


Ch^er 4 


CCMBUSTION IN A TOIN 
ANNULAR CHAMBER 

In the present chapter attention will be given to nonlinear 
instabilities in a thin annular combistion chand>er with uniform 
injection of liquid propellant (Fig. 3). The unsteady combustion 
response is characterized by a sin?)le tin»-delay mechanism which is 
the same as in a thin rectangular chandler investigated in Chapter 3. 

For this case the potential function $ depends only on 6 and t. 

Therefore, (2.36) can be written as: 

<l>,tt ■ * n((>^,^) + 2(j),g(|»,g^ - O(e^) * 0 

(4.1) 

This potential function predicts the behavior of the anplitudes 
of the various acoustic modes of the chamber. The velocity potential 
will be expanded in terms of products of the acoustic normal modes and 
time -dependent modal anplitudes as follcws: 

(}>C0,t) * Fj^(t)cosO + F2(t)cos23 + (4.2) 

This represents a solution involving standing waves. Traveling 
waves con also be treated with no essential difficulty, hut tliis matter 
is not pursued in the present work. 
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Substituting (4.2) into (4.1), ejqunding all nonlinear terms in 
a Fourier cosine-seties, and equating the coeffici^ts of cos 6 and cos 26 
individually to zero yields, 

fJ + Fj ♦ e(o((l - n)Fj^ ♦ nF^^) + 2(F^F2 ♦ F 2 FJ,) + kF) ♦ O(e^) - 0 

(4.3) 


F* ♦ 4 F 2 + e(o((l - n)F 2 + nF 2 q.) * * 4*cF2) ♦ O(e^) - 0 (4.4) 

These equaticais are roathanatically analogous to the equations 
obtained by Powell [10] in his invest igatiai of transverse instability 
in a cylindrical chamber. He obtained numerical solutions to his 
equations. In vdiat follows the methods of multiple sales will be used 
to deduce £q)proxiinate analytical solutions capable of explaining the 
results con^juted by Powell [10]. 

To begin the perturbation procedure let 

e - t (4.5) 

n = et (4.6) 

By using (4.5) and (4.6) eqiuitions (4.3) and (4.5) can be 
readily written as: 
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* ^^2 * ‘ *')**2*C * 


- F^Fj.O + O(e^) - 0 


(4.8) 


It will now be assuned that the first and secoid modal auplitudes 
can be e]q>ended in the fom 


h * ^iO ^ ^^il * "" 


(4.9) 


\diere i can be either 1 or 2. 

Substituting (4.9) into (4.7) and (4.8) yields 


^lO'CC ^10 ^ •'^lO’CC * ‘ 


+ n(Fj^Q.j.),^) + 2 (F^qF2q,^ + ^20^10,0 ■*■ ^ll’CC 


+ O(e^) - 0 


(4.10) 




- n(F20T),^) - FjqF^o.^ * . 4F2j) - O(c^) - 0 (4.11) 

Since Fj^q, F^j^, F-,q, and F2j^ do not depend on e each coefficient 
rf a power of c must be individually set equal to :ero. This yields 


F ^ F * 0 


(4.12) 




20 


0 


(4.13) 
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’ "^*’10*5 * 

+ 2(FjqF 20,^ ♦ *'20*’l0*5^ * *’ll*55 ^ ^’ll " 

^*^20’5n ‘‘‘ ^*^*^20*55 * ' *'^*’20*5 * *'^*^201^*5^ 

■ *’lO*’lO *5 “■ *’21*55 ^ ^*’21 " ° 


(4.14) 


(4.15) 


The general solutions of equations (4.12-13) may be vnritten as: 
Fio(^»n) “ Aj^Q(n)cosi5 + B^gSiniC (4.16) 

where i remains either 1 or 2. 

TTie coefficients A^g aiui B^g can be detennined by substituting 
(4.16) into (4.14) and (4.15), and equating the coefficients of cos5 
and sinC in (4.14) and cos25 and sin25 in (4.15) individually to zero to 
eliminate singular terms. This produces 


A^O + ^(a((l • n)A^g + n(.4j^gCOsT - B^gSinT) 

^10^20 ®10®20 * ‘'*'10^ “ ® 

Bio ♦ ^-(^((l - n)Bjg + n(AjgSinT B^gCOsT)) 
* ^10**20 ' "^!0^10 ‘ *^^10^ " ° 


(4.17) 


(4.18) 
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^0 ~ *'(A2qCOs 2T - v*”"» 

* ?t^0 - * **®2o’ ■ “ 

^20 * ^^C(2 - n)B2Q ♦ n(A2QSUi2T + B2 qCOs2T)) 

• Ko®ll) - ®^20> ■ 0 «-2») 

In obtaining (4.17-20) tl» identities 

cosCsin2C ■ 2^®^^ sin3C) 
cos^cos2C « yCcosC * cos30 
sinSsin2^ - -jCcosC - cos3f;) 

were used. 

It is com'enient to rewrite (4.17-20) in terms of amplitudes 
and phase angles defined as: 

Ai « C^cosa^ (4.21) 

“ C^sina^ (4.22) 

For additiaial convenience the final equations will be written 
in terms of C^, C,, and 

®i2 “ '^2 * ’“l 


( 4 . 23 ) 


Substituting (4.21*23) into (4.17-20) md rearranging 
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yields 



Cl + ^(oCj^(l - n(l - COST)) ♦ C^C2CosOj^2) “ 0 

(4.24) 


^1 * ■^(^2*^12 * - k) - 0 

(4.25) 


C 2 + i€(oC 2 (l - n(l - cos2T)) - ^^coso^2) “ ° 

(4.26) 


1 

^12 ^ ’ 2C2)sinaj^2 * on(sin2T - 2sinT) - 6k) « 0 

2 

(4.27) 

The linear stability analysis can be recovered from (4.24-27) by 
neglecting all nonlinear terras. Inspection of the resulting equatiem 
shows that the first mode is stable for 


1 - COST 

(4.28) 

while 

the second mode is stable for 



ITT cos'rr 

(4.29) 


It is desired to investigate situations for which the first mode 
linearly unstable vshile the second mode is linearly stal)le. That is, 
situations for which 


1 .... 1 
■c5iT-"- r” cos2T 


(4.30) 


r 
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For (4.30) to hold the inequality -y £?< -y oust hold. In such 
cases a limit cycle is possible. Its anplitiate can be obtaii^ from 
(4.24*27) by looking for solutions such that and C 2 are constants. 

It can be sem that such solutions exist only if 0^2 constoit (note 
that and thus 02 need not be constant). 

Solving the algebraic equations produced by assuming that 
C 2 » and 0^2 are coistants yields 

Cl - ^ 5 ^ / (1 - n(l - cos2T))(n(l - cosT) - 1) (4.31) 


C . o(n(l - COST) * 1) 
2 


(4.32) 


^ 6</o + n(2sinT - sin2T) 

tanOj^2 ■ 3 ♦ n(2cosT ♦ cos2t - S) 


(4.33) 


To determine whether or not the solutions of (4.24-27) actually 
approach the limit cycles discussed above, these equations were solved 
numerically by the Runge-Kutta method (no closed form solutiais being 
jqjparent). Because of the C 2 appearing in the denomination of (4,27) 
it is necessary to use a power- series solution for small times whenever 
the initial condition 0,(0) » 0 is used. The appropriate series 
expansions can be written as: 


Cl - 1 ♦ C^^t ♦ 0(t") 


(4.34) 


C 2 - C2^t ♦ 0(t") 


(4.35) 


I 
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Oh - * O(t^) 

(4.36) 

®12 * ®121* * 

Substituting (4.35-38) into (4.24-27) yields 

(4. 37) 

■ 1 - |oc(l - n(l - cosT))t 

(4.38) 


(4.39) 

Oi • -^(k - ^inT)t 

(4.40) 

♦ ^OT(2sinT - sin2T)t 

(4.41) 


In order to investigate some aspects of tte nonlinear stability 
pit^leni a ccmsiderable amcHint of minerical calculations were performed. 
Only the most inixirtant conclusicms and some typical results will be 
reported in the remainder of this ch^ter. These results are presented 
with the following objectives in mind: 

(a) The prediction of stable limit cycles for transverse mode 
instability. 

(b) The determination of the dependence of the amplitudes of 
tte limit cycles i:pon the caifcustion parameter n aid time-delay T wiiich 
are given by the steady state solution. 

The regioi of tte n - T plane in which limit cycles are expected 


is shown in Fig. 10. Defining 
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n, ^ (4.42) 

^ 1 - COST 

n . 1 (4.43) 

^ 1 - COS2T 

It can be seen that for (4.30) to be satisfied 

n^<n^i2 (4.44) 

Thus the region ABC (where B corresponds to n = ®) is the 
region of interest. Points A and C can be found fim n^^ = n 2 or 

cos2T - cosT = 0 (4. 45) 

In the region 0£T£2 it (4.45) has the roots T = -^ and 

To determine whether or not solutions to (4.24-27) actually 
approad'iOd the limit cycles predicted by (4.31-33) a variety of 
numerical solutions to (4.24-27) were obtained by using a fourth-order 
Runge-Kutta method. Sane typical results are shown in Figures 11-15. 

Figures 11 and 12 present results for T = 150°, n = 0.65, and 
< = 0. Substituting these values into (4.31-33) yields 


.^2 ~ 1.368 Rad. 

(4.46) 

= 3.768a 

(4.47) 

C 2 = 1.058O 

(4.48) 


for the limit -cycle values. 
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It can be sem that tlw numerical results do approach the 
limit -cycle values and that t!»se are in agreanent with (4.46-48). It 
was established that other sets of initial ccmditions lead to the same 
limit cycles. The results of the con;>utations used to establish this 
fact are not shown. 


Figure 13 shows some typical results for T * it 
this case it is possible to siinjlify (4.24-27) to 

and K ■ 0. In 

• 1 

Cl + |e(oCj^(l - 2n) + Cj^C7Cosaj^2^ “ ® 

(4.49) 

• 1 

+ |eCj^C2sinaj2 “ 0 

(4.50) 

* 1 12 

C., + ^(oC, - ^ Cj^cosaj2^ * ^ 

(4.51) 

“12 ■" ■ 2C2)sina^2 " ^ 

(4.52) 

For the initial conditions aj^(O) » 0^,(0) * 0, 
are satisfied by 

(4.50) and (4.52) 

“l ' “12 ' “ 

(4.53) 

Substituting; (4.53) into (4.49) and (4.51) yields 

+ ^(o(l - 2n)C^ + CjC2) = 0 

(4.54) 


- kh = 0 

it Atf H X 



(4.55) 
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Hius and ^ determined in close- form and and C 2 

must be determined numerically from (4.S4) and (4.S5). For these 
conditions (4.31-33) simplify to 


' 2o Vln - 1 

(4.56) 

' o(2n - 1) 

(4.57) 

“12 “ ® 

(4.58) 


It can be seen that (4.56-57) are consistent vdth (4.54-55). 

The numerical results shown in Figure 13 are for n ■ 1. It can 
be seen that the correct linut-c>'cle values are approached for large t. 
No difficulty was encomtered in getting the numerical solutions to 
approach the limit-cycle values given by (4.54-55) for any value of n 
tried. This is in ccMitrast to the situaticms for T n. This will be 
discussed below. 

For values of T other than ISO® there appeared to be an upper 
limit on the value of n for which the nunerical solutions would approadi 
the limit -cycle values given by (^4.31-53). (The upper limit appeared 
to be approximately = O.'i’S, but more work would be required to imike 
this value more definitive.) For n>n^^ numerical calculations indicate 
that Cj - :uid ‘ as t ‘ Inspection of the nmnerical results 
showed that Cj and C, would cene close to the limit -cycle values for 
intennediate values of time, but then diwrgc from these values because 
ivould pass the limit-c>’cle values and ccMitinue growing. Thus it is 
suspected that n^^ “ ^ for T » 180° and k; * 0 because Uj, (t) - 0 for 


41 


these coiditions. Thus is already at its limit-qrcle value and 
r«nains there for all t. It is suspected that the limit cycle is 
iBistable for n>n^. A linear stability analysis could be carried out in 
an atteiqpt to verify this conjecture, but this was not draie in the 
present work. 

To make sure that the lack of approach to the limit cycle 
discussed above was not due to errors in the mmerical procedure, the 
following method was used. A value of selected, and (for given 

values of T, a, and n>n^) the correspond value of k was found fran 
(4.33). Then the appropriate values of Cj^ and C 2 were found from 
(4.31-32). Then a nimerical solutiai was found by using these vali;«s 
of Cj^, C 2 , and Oj ^2 ^ initial conditions. It was found that in 
these cas^s the variables remained ccmstant at their initial values. 
Results of tw t)’pical calculaticms of this type are shown in Figures 14 
and 15. It was also observed that the slighest deviation of the initial 
values from those confuted by the procedure described above lead to 
divergence of the variables fran their limit-cycle values. This 
behavior provides further support for the conjecture that the limit 
cycles are unstable for n>n^^. 

It should be noted that for k = 0(1) the frequency of the second 
acoustic mode is very nearly twice the frequency' of the first acoustic 
mode. The perturbation method makes it clear that the existence of a 
limit cycle depends critically on this fact. This is not evident from 
the numerical solution found by Powell [10]. 


Chapter 5 


CONCLUSION 

In this report the two-variable i«rturbation method was used to 
solve the modal -ai^litude equaticms arising from the application of the 
Galerkin method to a nonlinear wave equation describing pressure-sensitive 
combustion- instability problms in liquid-fuel rocket motors. Two 
idealized problems were discussed. The first was that of transverse 
synmetric motion in a narrow rectangular cranhistion chamber. This problem 
exhibits self-coupled instability. Tl;e second was that of transverse 
motion in a narrow annular comtaistion chamber. This problan exhibits mode- 
coupled instability. In this way an application of the method to each type 
of pressure- sensitive instability problcan was presented. 

For the case of the i octangular chamber both the linear and nonlinear 
stability boundaries and limit-cycle amplitudes were obtained in closed form. 

It was also possible to determine the evolution of the initial disturbance 
in closed form after making a slight approximation. When this approximation 
was removed a numerical solution was required to determine this evolution. 

For the case of the annular combustion chamber the linear stability 
boundaries and the limit-cycle amplitudes were foimd in closed form. Tlie 
evolution of the initial (Jisturb;incc was determined luuncrically. 

In both cases the results exhibited the s;ime behavior as tlvit observed 
by Powell flO] who used direct numerical solution of the modal -amplitude 
c\]uations to ;uialyze similar problems. In fact it w'as found tlvit the closed- 
form solutions obtainexl herein could be used to explain several of the 
pariunetric trends obsen^ed in [10]. Since the functions determined numerically 
in the present work arc much simpler than those compH-itcd in [10] it is thought 
that present procedure offers opportunities for significant savings of 
computer time. 
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0(1 the basis of the nianerical results the following nonlinear 
mechanisms were found to be important in determining the nonlinear 
stability characteristics of the syst«n: (a) the transfer of energy 

between inodes and (b) the self coipling of a mode with itself. 
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Figure 4. Stability Boundaries for Rectangular Chamber (exact theoiy) 










Figure 7. Modal Airplitudes Versus Time for Rectangular Chairijer 
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Anplitudes and Phase-Angle Difference Versus Tijne for Annular Chamber 









Figure 13. Modal Amplitudes and Phase-Angle Difference Versus Time for Annular Qiamber 
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Nbdal Airplitudes and Phase- Angle Difference Versus Time for Annular 



Figure 15. Nkxlal Amplitudes and Phase-Angle Difference Versus 






